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We present results for the local density of states in the S and N layers of a SN multilayer, and 
the supercurrent, based on a Green's function formalism, as an extension of previous calculations 
on NS, SNS and SNSNS systems. The gap function is determined selfconsistently. Our systems are 
chosen to have a finite transverse width. We focus on phenomena which occur at so-called critical 
transverse widths, at which a new transverse mode is starting to contribute. It appears, that for an 
arbitrary width the Andreev approximation (AA), which takes into account only Andreev reflection 
at the SN interfaces, works well. We show that at a critical width the AA breaks down. An exact 
treatment is required, which considers also ordinary reflections. In addition, we study the influence 
of an interface barrier on the coupling between the S-layers. 

I. INTRODUCTION 

Starting about two decades ago, the interest of developing devices at a very small scale gave rise to a new branch in 
physics, the mesoscopic physics. Both theoretically and experimentally, many interesting phenomena were discovered 
which occur at this scale which lies essentially in submicron ranges. 

Many samples are built up out of superconducting (S) and normal metallic (N) components in which necessarily 
SN interfaces and possibly point contacts occur. It is why a lot of theoretical work was devoted to studying different 
SN configurations. 

The first experimental investigations, using tunneling spectroscopy measurements, revealed the fact that the density 
of states in a normal metal connected to a superconductor is modified 1 . McMillan provided a simple tunneling model 
for the proximity effect at SN interfaces, which allows for a solution of the Gor'kov equations. Ishii 2 and Furusaki 3 ' 4 
extendend his work to include Andreev reflections. 

Recently 5 a powerful Green's function formalism was published which unified earlier formulations 1-4 and improved 
upon them. First applications were made for NS, SNS, and SNSNS systems 5-7 . An important feature of these 
calculations was, that the systems were chosen to have a finite transverse width and they were focussed on particular 
phenomena which occur at specially chosen transverse widths. Up to now only treatments are known for an infinite 
width 8-10 , or if a finite width was considered it was done in a global way, in terms of the number of allowed transverse 
modes 4 ' 11 ' 12 . 

The aim of the present paper is to show applications to different SN multilayer structures, by which we understand 
a periodic sequence of S and N layers, extended in the x-direction, see Fig. 1. A set of observables is calculated, such 
as the local density of states (LDOS) in the S and N layers of a multilayer, and the supercurrent in the multilayer. 

In Section II, we present the theory, applied to SN multilayers. In many situations, it is enough to work within the 
Andreev approximation (AA), which reduces to taking into account only Andreev reflections at the S/N interfaces. 
In Section III we will show results derived within AA. Exact calculations, which include ordinary reflections, are 
discussed in Section IV, in relation with the so-called critical transverse widths, at which a new transverse mode is 
starting to contribute. In Section V we present results for the supercurrent in the SN multilayers. In Section VI we 
study the consequences of using a selfconsistently calculated gap. Finally, to complete the picture, in Section VII we 
take into consideration an interface potential, by this modelling a Schottky barrier. We apply a simple (5-function 
barrier located at the SN interface, introduced by Blonder et al. 13 . 

II. CLEAN SN STRUCTURES 

The purpose of this section is to summarize first the basic ingredients of the general theory presented in Ref. 5 
which we need in calculating the LDOS p{x, E) for energies E of the order of magnitude of the gap energy and the 
supercurrent /. After that we show how the general theory is elaborated for applications to SN multilayers. 
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We study a periodic SN multilayer, which extends in the x- direction, as depicted in Fig. 1. In the transverse 
directions y and z, the system has finite size, L y — L z — L t . We apply a Kronig-Penney superlattice model, depicted 
in Fig. 2, which means that the pair potential is 

A(x + d s + d N ) = A(x)e** 

, _ J A ifxeS layer , . 

[X) ~ \ ifxeNlayer, [ ' 

where ds and are the thicknesses of the S and N layers respectively, and <fi is the phase difference between two 
neighbouring superconducting layers. 

Both quantities we want to calculate, can be expressed in terms of the Green's function in the following way 

2 1 x 

p(x,E) = lim } lmGu(x,x;ky,k z ,E + iS). (2) 

7T L y L z S^O 



2ie ^L z £ X H - kT^Gnix, x'; k y , fc„ 



(3) 



where the Green's function Gu is the upper left element of the matrix Green's function. 

The Green's function can be expressed in terms of a kind of wave function, which is the solution of the one- 
dimensional differential equation 

/ • d 2 ,2 

iu n + — + k Fx -A , 

dx 2 d2 |*(x) = 0. (4) 

where k 2 Fx = fi — k 2 — k 2 . Note that the Bogoliubov equations arise from Eq. (4) by substituting iui n — > E. The 
solution of Eq. (4) for the spatially dependent A of the multilayer (1) is built up by starting with the solution for a 
homogeneous superconductor, having a constant Ag. The wave function can then be written as 

u -* e -i<t>s/2 J 



with tig = \JiLo n + iafls, i^s — \J + 0+) 2 — A|, and the Matsubara frequencies u) n = rinksT, n taking odd 
integer values only. The four standard solutions are labeled with the sign indices a and v, that can both equal ±1. 
The index a refers to the type of the propagating particle (electron-like for a = + and hole- like for a = — ) and the 
index v indicates the direction of propagation. 

In order to express the Green's function in terms of ^^(x), a conjugate wave function is needed, namely, 

*g"(x)=^ uge-** 8 / 2 Ug "e**"/ 2 ^ e iavklx (6) 

(which is not the hermitian conjugate), in which v has now to be explained as minus the direction of propagation. 
With the use of these wave functions, we can express the Green's function for a homogeneous superconductor as 

g s (x,x')=j2d s ^(^rnx% (7) 

a 

with /u = sgn(x - x') and d a s = - in \ k a - 

A. A single interface and more interfaces 

For a single interface, situated at the position Xj , the general form of the Green's function is 
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G v j v ij{x,x') = g vj (x,x')5 vv > +Y,dl j di,^{x)tlfjf^\x'), (8) 

aa' 

where the subscript (uj) refers to the part of the system that is at the v side of the interface at position Xj. The first 
term accounts for the possible ways of propagating from x' to x without being scattered at the interface, while the 
second term describes the propagation via the interface. 

The scattering matrix t°J v Yp is found by applying boundary conditions at x = xj . We require the continuity of the 
Green's function and its derivative. Thus, the equation obeyed by the i-matrix is 

Y.^i^M^ii't = -srfr'fo)- ( 9 ) 

If wc consider a system with an arbitrary number of interfaces, with position coordinates Xj < Xj + i, then the 

scattering of the quasiparticles is described by the scattering matrices T^^.f . The general form of the Green's 
function is 

Gi/jif'j' x ) — Q v j{x^ x ) [Sw'Sjj' ~\~ &—w'&j-\-v,y\ 

+ E E dZjti^ixKfjF' (io) 

aa' [ifj,' 

Again, imposing the boundary conditions, we obtain a Lippmann-Schwinger equation which allows to calculate the 
T-matrices by means of the single interface i-matrices 

T vjJj> =t lji^j [V"'^J' + &-n'v'?>3+li',j'] + E d l"j T v"ju ! j ,V ^ • (I- 1 ) 

<t"v" 

This equation expresses the idea of multiple scattering, since the matrix T^Jy contains all possible processes that 
yield the correct final state. The first term in Eq. (11) accounts for the possibility that the particle is scattered once. 
The second term collects the processes in which the particle is scattered once due to i^/JJ and an arbitrary number 
of other times due to T°„j v !j, v p . 



B. Periodic SN multilayers 

Till now we just summarized the description given previously 5 . We will now focus on an infinite periodic SN 
multilayer to which the theory has not been applied yet. For an infinite multilayer, it is always possible to refer to 
any layer of the system by referring to an even-numbered (or an odd- numbered interface only). This can lead to a 
further simplification of the Lippmann-Schwinger equation. We choose to refer to any part of the system by referring 
the even interfaces. However, for the i-matrices, both even and odd interface indices need to be used. 

In order to rewrite Eq. (11) for the present purpose, we define the following matrices 

rpaa rpaa , — fj,/j, 

1 rpaa , — rpaa [ifi 



-M^'j' -ni,-y'i' 



Aj = I |j C -M,j+M.f*'>j+M a -M'>j-2°M-V- I (14) 

, u t-nd+w'j+^-n'jO-w' (15) 

,aa ja n / 

t -M i i'i Vi U / 

( 4f>°' -W' r ia' r s \ 

C 3 ;=[ Z - l i,3+y.,y', 3 + l i a - l i', 3 +2 t>+ y'+ (16) 
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All elements of these 2x2 matrices can themselves be regarded as 4 x 4 matrices, with the indices (a, fi) and (a 1 , ji') 
labeling the rows and the columns, respectively. That makes TV an 8 x 8 matrix that satisfies 



A 3 ' T J-2,j' + ( B J - 1) ' TV + C j ■ T 3+2,j> + Vj5jj> - 0. 

In terms of these matrices, the system Green's function (10) gets the following form 

G v j v iji{x,x') — G V j{x,x') [5 vv i5jji + 5- vv i5j +v ji] 



(17) 



(TV - ^MVJt' (■■■A j ,B j ,C j ,A j+1 ,B j+1 ,C j+1 , (x% 



(18) 



All possible scattering processes are incorporated in Eq. (17), which expresses the content of the Lippmann- 
Schwingcr equation. To illustrate the different processes accounted by the Aj, Bj, and Cj matrices, we show a 
schema in Fig. 3. 

We now turn to the periodic system. For the moment, we assume that the phase of the pair potential is the same 
and equal to zero in all the S-layers. The periodicity allows us to simplify the problem and to rewrite equation (17). 
First we perform the following transformations 



''vjv'j — e 'ViV j e ' 



(19) 



rpaa'uu' _ lavK-^jXj rpaa vv 10 f n^^^j 



(20) 



The scattering matrices with hats no longer depend on the interface positions Xj, although they still refer to the 
interface number, through the labels vj. 

The matrices TV , D ^ , Aj , Bj , Cj , become 



T jj' = D j%' + 



rpCTCr' , — flfl' rp<T<j' , — fifi f 



D 



(7° 



V 







d 



7T c 



U 1 a - i i'°ii-°\x'- e M 



ftrcr' \-nn' 
l i l- i l' a -i 




B 







faa'-nn' c 



i& f k a _ ,a_ 



t -nn' d n' e M 

n faa' -\i\i' , a ' r r ia'k°',a_ , 





(21) 

(22) 

(23) 
(24) 
(25) 



After these transformations, and as a consequence of the periodicity, the Green's function becomes dependent on 
the relative coordinates only, and the j-independent layer thicknesses a M = [i(xj+n — Xj) and a-^ = fj,(xj — Xj-^) can 
be defined. In terms of these matrices, Eq. (17) now reads as 

A • Tj- 2J , + (B - 1) ■ T jr + C ■ T j+2J , + t)6 jr = 0. (26) 
This is a kind of discretized version of the original Green's function equation. The general solution is 

- *££-?') • f 0' ( 2? ) 

where X_, X + and T arc implicitly given by the following set of equations 
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C = 

2 



A • X 2 + (B - 1) • X_ 
A + (B — 1) ■ X + + C • X+ = 
[A • X + (B - 1) + C • X+] • f o + D = 0. 



(28) 
(29) 
(30) 



Note that (28) and (29) are quadratic equations. By that, for a periodic system, we can rewrite the expression of the 
Green's function (18) in a simpler form 



G vjvlj .{x, x 1 ) = Q vj {x, x') [cW%' + 5- vv '5 j+Vtj ,\ + d vj d °''j'Kj( x - x i) 

(7(7 ! f_L/_L f 

(t jy - D,%0^'(A,B,C)*^'(x' - x,,), 



(31) 



The problem of calculating the LDOS or the supcrcurrent, reduces to solving the system of quadratic matrix equations 
(28) to (30). 

After some manipulations, one can manage to reduce the problem to solving a 2 x 2 matrix equation, which is 
equivalent to solving a system of 8 simultaneous equations with real coefficients. In the Appendix, we show this more 
explicitly. 

The formalism described up to now can be extended to the situation in which the pair potential has a phase, which 
allows for currents in the multilayer. 

By convention, we assume that the N-layer has also a phase, equal to the phase of one of the two adjacent S-layers. 
This means that the phase over a bilayer is constant and it makes a jump of <f> at the edges between two bilayers. 

Suppose the interface between two bilayers in chosen at j odd, then for even j the i-matrices obey the equation 



(32) 



while for odd interfaces, a <f> dependence remains and we are left with 



(j odd), 



(33) 



where 



U{<f>) 



e^/ 2 
e-^/ 2 
e^/ 2 
er^l 2 



(34) 



III. LOCAL DENSITY OF STATES IN SN MULTILAYERS 



First we apply the theory in calculating the LDOS in the middle of one of the S or N layers. For most of the systems 
which will be described, the transverse width L t is fixed to 13 Bohr, the chemical potential fi = 0.5 Ryd, and the 
LDOS is normalized to the spatially constant LDOS of the bulk N material. The coupling potential V is calculated 
using the BCS formula 

T c = 1.13uj D e- 1/NMV , (35) 

where N(fj.) = ^/JI/4ir. For Al with T c = 1.2 K, /j, = 0.5 Ryd, and oj d = 375 K, we find V = 9.516 Ryd. For the 
pair potential A we choose a value of 0.0001 Ryd. Given the BCS relation -^p- = 1.77, the pair function should 
be somewhat larger, A = 0.00018 Ryd. However, we should also keep in mind the reduction of the pair function 
due to the finite size of the system. This we will discuss in Section VI, in which we will determine the gap function 
selfconsistently. Furthermore, many of the results we will show are not quite sensitive to the precise choice for A. 

It appears that the LDOS curves for SN multilayers look rather complicated. As a preparation to understand them, 
we first look at simpler systems and we postpone the treatment of SN multilayers to subsections IIIB and IIIC. In 
the coming subsection, we will follow the development of the LDOS for a bulk system to the LDOS for systems with 
a few interfaces. 
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A. From bulk superconductor to SNS system 



In this subsection we will first follow the development of the LDOS for a bulk bar-shaped superconductor to the 
LDOS of a SN multilayer with ds » (£,c(/v) arid of multilayers with d$ » £, in which becomes comparable 
to ds- For the present clean systems the BCS coherence length £ w 4000 Bohr. Looking at Fig. 2 it is clear that 
a multilayer with thick superconducting layers, such that ds >> £, comes close to a SNS system, particularly for 
E < A. 

In Fig. 4 we show the LDOS inside the S layer of a SN multilayer and the DOS of bulk S material. One clearly sees 
the singularity in the LDOS at E = A. The non-zero DOS of the S material for E < A is due to the small imaginary 
part iS added to the energy, E + iS. In all calculations we used S ~ 0.02A. Here dN = 1000 Bohr and ds = 50000 
Bohr, so that the presence of the N layer is just a small perturbation from a bar shaped superconductor. 

First we concentrate on the development of the LDOS for E < A. In Fig. 5 we show for the same system both the 
LDOS inside the S and the N layer. Due to the very small N layer thickness, there is just one Andreev bound state 
in the N layer LDOS close to the gap value of the energy, which is broadened by S to a peak. 

Figs. 6, 7, and 8 show what happens to the LDOS of the SN multilayers if we increase dN to dN = 2000, 4000, 
and 10000 Bohr respectively. The pictures look more and more complex as we increase djv In the N layer LDOS we 
notice the appearance of more Andreev bound states at lower and lower energies. The singularity in the S-layer DOS 
lowers, certainly due to reduced interaction of the neighbouring S-layers. 

The oscillations in the LDOS in all figures for E > A are a periodic-multilayer effect. Before discussing this 
band structure effect, we will first consider multilayer systems with decreasing ds, by that making more explicit the 
multilayer character of the system. 



B. From SNS system to SN multilayer 

Let us pick up the N layer LDOS from Fig. 8 and put it together with the LDOS of a SNS system, whose 
dN = 10000 Bohr. This is what we show in Fig. 9. The good similarity is due to the large ds- This SN multilayer is 
just a perturbation of a SNS system, as we can notice from the small oscillations at energies E which are larger than 
the gap. Way below E = A the Andreev bound states curves coincide. Just below E = A the clear peak in the SNS 
system curve is smeared out in the multilayer curve, due to tunneling interaction between the N-layers. 

The similarity between SNS systems and SN multilayers reduces with decreasing ds- Multilayer features start to 
appear gradually in the LDOS, as we can see in Figs. 10 and 11. 

Due to the increased tunnelling, the discrete states start to form bands, while at E larger than the gap, the 
oscillations are more pronounced and follow a periodicity, according to the dispersion relations 8 



ui. i. ±j w ,VE^+Wds jEd N 
cos[{k x - k Fx ){d s + d N )\ = cosh — c oshl ^— - + — 

VWT^ s jEdN + i± 

2k Fx I 2k Fx 2 



cos[{k x + k Fx )(d s + d N )] = cosh ^^ ~|~ — ^ S cosh 

lk Fx 




u ^WTWdl^ ( Ed N i<j> 

where k Fx is given by 



+sinh v ^ o 7 1 " 1 " fa sinh - ^ , (36) 

2k Fx \2k Fx 2 J' K ' 




In Fig. 12 we illustrate these dispersion relations for a multilayer with ds — dN — 10000 Bohr and two different 
choices for the phase difference, (f> — and </> = ir. To make it more clear, we considered here only the (l,2)-mode 
contribution to the LDOS. We notice the change in the succession of gaps with the phase 4>. 

Using the definition of k Fx , we notice that the number of allowed modes (n y , n z ) is limited by the condition k Fx > 0. 
Besides, since k Fx is different for each mode, the periodicity with which there is a solution for k x in the equations 
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(36) is also different. For a higher mode (n yi n z ), kFx is smaller and the frequency with which peaks and gaps in the 
LDOS are alternating increases. We can see this in Fig. 13, where we compare the contributions to the LDOS coming 
from different modes, at phase 4> — 0. As in Fig. 12, the system has ds = (In — 10000 Bohr. One can easily check 
out that by just adding the four modes' contributions one obtains the total LDOS shown in Fig. 14, which we are 
going to discuss in the following subsection. 

C. From infinite transverse size of a SN multilayer to a finite one 

Now we are prepared to investigate real periodic-multilayer effects. In addition to the (p — result shown in Fig. 
14 we also show the 4> = ir result for the same system in Fig. 15. The dispersion relations (36) are nicely illustrated 
in the succession of bands and gaps. This system was studied also by Tanaka and Tsukada 8 . In their Fig. 2a, they 
describe a SN multilayer which has d,N = ds = 5000 A and <j> = 0, and is infinite in the transverse direction. For 
E < A the pictures look quite similar, but above the gap the LDOS in Figs. 14 and 15 is much less smooth than in 
Fig. 2 of Tanaka and Tsukada. However, in Fig. 16 we show the LDOS for a system with a larger transverse width 
(L t — 130 Bohr), and indeed, the behaviour for E > A has become much smoother. So we conclude that differences 
between our results and the results of Tanaka and Tsukada come from their use of an infinite transverse width. 



IV. EXACT CALCULATIONS AT CRITICAL WIDTHS 

In many situations, the Andrcev approximation gives results with an error which is estimated to be less than 0.1%. 
The systems which we discussed up to now satisfy the conditions for which the Andreev approximation is very good. 
In this section, we will deal with situations in which exact calculations are necessary. 

In Fig. 17 we show the not-normalized LDOS at E = 5A of a homogeneous bar as a function of the transverse width 
L t . At E = 5A, the LDOS of multilayers with the same L t approaches a constant value, given by the homogeneous 
bar, as we can see in Figs. 4 - 16. We notice in Fig. 17 that at certain values of L t , the LDOS has steps, followed 
by a fast and smooth decrease. At these widths, where the condition k Fx > reaches the equality, an extra mode is 
allowed in addition to the previous ones. The new mode has a large contribution to the LDOS, explaining the step. 
Apparently, at smaller L t 's the steps are higher, which means that the effect of adding a new mode is larger. These 
values of L t are called critical widths, and we will denote them by Lf. 

Close to the critical widths, for < k Fx < A, the AA is not good anymore. The highest modes contribute most to 
the LDOS, as the steps in Fig. 17 suggest us. Besides, since kFx is very small in the dispersion relations (36), these 
modes will give rise to much more states than the lower modes. 

As an illustration, in Fig. 18 we show the absolute value of the LDOS for a SN multilayer, whose ds = = 10000 
Bohr, (like we showed in Fig. 14), but this time at a width L t = 12.566371 Bohr. This value of the transverse 
width, corresponding to the equality k Fx — 0.3A for the highest, (2,2) mode, lies very close to the critical width 
corresponding to k 2 Fx — and can be implemented without getting numerical problems. At this width, we calculated 
the LDOS at two different positions x with respect to the SN interface, inside the N layer. Although the solid and 
dashed curves have peaks at the same energies (the dispersion relations do not change with x), their magnitude goes 
up or down, depending very much on x. This is not the case within the Andreev approximation, represented for 
comparison with a dotted line, where peaks of the same height are situated on the energy axis at equal distance from 
each other. This difference between the exact and the approximate results can be explained if we make use of the 
definition of the LDOS, 

LDOS(r) = J2\*n(r)\ 2 S(E-E n ). (38) 

n 

In the AA the Andreev states with the electron moving to the right and to the left are uncoupled and are degenerate. 
They can be represented by plane waves, having a r-independent absolute value. In the exact treatment the corre- 
sponding travelling waves are coupled and they are split into two standing waves, an odd (sinus) and even (cosinus) 
function. This leads to weighting factors in the expression of the LDOS which are different and position-dependent. 
So, the lifted degeneracy in the exact calculation explains the position dependence of the LDOS illustrated in Fig. 18. 

In order to compare multilayer results with publicated SNS results 5,7 , we calculated the LDOS inside the N and S 
layers of a multilayer with djy — ds — 4000 Bohr. This is shown in Fig. 19 together with the comparing SNS result. 
The corresponding results derived in the Andreev approximation are shown in Fig. 20. We restricted the calculations 
to the highest mode's (2,2) contribution to the LDOS, at L t = 12.5676 Bohr. This transverse width corresponds to 
k Fx — A for the mode (2,2). In both Figures, at energies E < A, we do not notice any difference for the N-layer. The 



7 



multilayer features appear only above the gap. Besides, inside the S-layer there is no contribution from the highest 
mode, as the corresponding states have such a small momentum k x that, for E < A, they are localized inside the N 
layer. 

The features shown in this section are directly related to a fine-tuning of the transverse width. In this respect, 
these results are new compared to those reported by Tanaka and Tsukada 8 , who have considered an infinite transverse 
width only. 

V. CALCULATION OF THE SUPERCURRENT 

This section is devoted to the supercurrent in a SN multilayer. In 1962 Josephson predicted that a supercurrent 
can be present in a SIS junction (Josephson junction) in the absence of an external voltage (dc Josephson effect). 
This current appears provided there is a difference (f> in the phase of the pair potential between the two S layers of 
the Josephson junction. 

I = IrnaxSm(j). (39) 

Further, if an external potential is applied to the junction, then 

f-W/l (40, 

In other words, an external potential gives rise to an alternating supercurrent of frequency / = 2eV/h (so called ac 
Josephson effect). The quantum energy hf equals the energy of a Cooper pair transferred across the junction. It 
appears that the Josephson effect is present also in SNS junctions. We will investigate it for SN multilayers. 

Using equation (3), we calculated the supercurrent I through a SN multilayer with ds — d^ — 10000 Bohr as a 
function of the phase difference <f> between two consecutive S layers. Fig. 21 gives the supercurrent normalized to the 
basic supercurrent unit Io — eA/h for different choices of the transverse width L t of the multilayer. The 4> dependence 
of I is basically similar to a sin0 dependence, as for a Josephson junction, in that it is periodic in 2tt. We notice that 
the supercurrent increases in magnitude with the transverse width, which makes sense, given the fact that the larger 
the width L t , the more modes contribute to the current. However, a small deviation from this monotonic behaviour 
in the dependence of the supercurrent on the transverse widths is noticed at a larger phase, <j> w 27r/3. In our Fig. 
21 we see this behaviour between the critical widths Lf = 14.0492, when the mode (3,1) starts to contribute and 
Lf = 16.0186, when the mode (3,2) appears. Around <j> = 2w/3, the curve corresponding to L t =15 Bohr lies slightly 
higher than the curve for L t —16 Bohr. This can be interpreted as being due to a destructive interference between 
the electronic contributions to the current at larger phase. This results in a small deviation from the symmetry of 
the sin-function dependence of the supercurrent as a function of phase. For an infinite transverse width, Tanaka and 
Tsukada show a similar dependence in Fig. 4 of their paper 8 . 

The way in which the transverse width influences the maximum of the supercurrent I max is shown in Fig. 22. 
The monotonic increase of the supercurrent exhibits steps at each critical width. This is not surprising, since at a 
critical width new modes start to contribute. At the onset of this contribution, the kinetic energy of the new modes 
£,2 _ ( "»'2°" 7r ) 2 — ( " z, ^ ax7r ) 2 is very small and so is their contribution to the supercurrent. But with the increase of 
L t , the supercurrent reaches a constant regime, till the next Lf. 

Now we fix the transverse width at L t = 13 Bohr and we change the layer thicknesses ds and djy. The results are 
shown in Fig. 23. If ds S> £, in which the coherence length £ sa 4000 Bohr, the SN multilayer compares well to a SNS 
system, for which 4> max ss 0.8-7T, as we will see below in discussing Fig. 30. At smaller values of ds the phase <j> m ax 
at which the current has a maximum shifts gradually towards lower values. Further, if the ratio between ds and djv 
is constant, the systems have approximativcly the same maximum supercurrent. However, it should be noticed that 
all systems have the same A. This picture of constant I max changes if the gap function is calculated selfconsistently, 
as will be shown in Section VI. With the decrease of d^ with respect to ds, the current increases due to a better 
coupling between the S layers. Modifying the value of djv doesn't have consequences on (j)rnax- This can be noticed if 
we compare the curves corresponding to d N = d s = 4000 Bohr and d N = d s /2 = 2000 Bohr. 

VI. SELFCONSISTENT CALCULATION OF THE GAP 

The formalism described in Section II can be applied to a selfconsistent calculation of the gap function A. The 
method, which is extensively described in Refs. 5,7 , is based on the selfconsistency condition 
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A(x) = -k B TV{x)^ ^G l2 {x,x;k y ,k z ,iLo n ). (41) 

* k y ,k z >0 n 

The summation in equation (41) is divergent. In order to render the sum convergent, a cut-off of the summation 
over the Matsubara frequencies is introduced, as in the following expression 

A(x) = -k B TV(x)jj J2 G 12 (x,x;k y ,k z ,iuj n ), (42) 

* k y ,k z >0 \u n \=wk B T 

where u>d is the Debye frequency. We limit the Matsubara frequency to lob> = n max irkBT ', where n max = [Od/7tT] 
and 9d = Wo/feg the Debye temperature. As we notice, at large temperatures, n max becomes small, while du) n — 
cj n — w n -i is large. This gives rise to big, unphysical oscillations of the order parameter A with temperature T, close 
to T c . We get rid of these unwanted oscillations by taking an integration, rather than a summation over u) n . Results 
for a bar of transverse widths L t — 30 and 100 are shown in Fig. 24 and Fig. 25 respectively, at which we will come 
back later in this section. 

In addition to the integration over the Matsubara frequencies, we also investigated another cut-off method, which 
avoids the gap oscillations at large temperatures. An alternative way to render the summation (41) convergent, is to 
impose the cut-off on the momenta k, instead of on the Matsubara frequencies. For a homogeneous S bar this reads 

k B TV(x) 1 rt*+ k B®D 

A{x) = Y2 / dk x V 'Y\G 12 {k x ,k y ,k z ,iu n ). (43) 

8^ L t J k 2 =fi—k B 6>D ~ 

More explicitly, equation (43) can be written 



[) ~ Stt L\ J k , = ^ kBeD dkx fa 2^ (i^ja — (fj, — k% — k y — klf - A» ■ 



We further perform the summation over the Matsubara frequencies, obtaining 

(45) 



V(x)A 1 [» +kB6 ° „ ^ tanh^ 



a(x) = — / dkx J2 

OK L t J k 2 = ^_ kB@D k - 



2E k 



where E k = (p — k x — ky — k^) 2 + A 2 . For bulk superconductors, both ways of rendering the integral convergent 
lead to the same result. However, in the case of a homogeneous superconducting bar, we have summation over the 
transverse momenta k y and k z instead of an integration. This affects dramatically the results when k y and k z are 
large, particularly at small L t . We can see this in the dependence of the gap A on the transverse width L t , shown in 
Fig. 26. The calculation is done at T = 0. The dotted line comes from a calculation with a cut-off on the Matsubara 
frequencies. Apparently, the latter cut-off method leads to a much more stable result than the method in which 
the momenta are cut off. The solid line exhibits unphysical oscillations in the gap indeed. In addition we show a 
dependence A(T) in Fig. 27. This curve exhibits a largely reduced superconductivity compared to the upper curve in 
Fig. 24 obtained by the other cut-off method. We conclude that cutting off the Matsubara frequencies leads to much 
more reliable results. Otadoy et at 5 - 7 used this method to calculate the selfconsistent gap for systems such as SNS 
and SNSNS systems. Here, we extend the application to SN multilayers. We show results for two transverse widths. 

In Fig. 24 the temperature dependence of the gap is shown for L t =30 Bohr. The solid, dotted, and dashed curves 
represent A(T) for a homogeneous S bar, a SN multilayer with ds — d^ — 10000 Bohr, and a SN multilayer with 
ds = dpf = 4000 Bohr respectively. As expected, the selfconsistent gap decreases with the periodicity ds + g?at of 
the multilayer. Indeed, for a smaller ds, the contribution to the averaged gap over the layer comes mostly from the 
regions close to the NS interface, where the suppression of the gap is most effective. Similarly, in Fig. 25 we show 
results for systems with L t =100 Bohr. Again, one clearly sees the suppression of superconductivity by reducing the 
transverse width. 

In addition to the results derived in Section V, in the present stage we can look to the temperature dependence of 
the supercurrent. First we show in Fig. 28 the phase dependence of the supercurrent for multilayers with L t = 30, 
at different temperatures. For a given layer thickness, the peaks at different temperature occur at the same phase. 
The multilayer with a smaller periodicity ds + d^ has the corresponding maximum at a lower 0, as we discussed 
in the previous section. We notice the suppression of the current with the increase of T. For the multilayer with 
ds = gJat=10000 and L t =30, we show the current-temperature dependence I(T), in Fig. 29. The temperature at 
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which the superrcurent becomes zero coincides with the critical temperature at which the corresponding gap function 
is zero, see Fig. 24. 

Finally, we compare the phase dependence of the supercurrent for SN multilayers with corresponding results for 
the SNS system. In Fig. 30 we show results at temperature T = 0.4 K. Obviously, since the gap function in the 
SNS systems is equal to the homogeneous bar gap, the supercurrent is larger than in the SN multilayers. However, 
with the increase of the N-layer thickness, the supercurrent in the SNS system decreases, since to a thicker N-layer 
corresponds a weaker coupling between the two half-infinite S-layers. On the contrary, for the SN multilayer results 
shown, the S-layer thickness increases as well when the N-layer thickness is increased, and the supercurrent increases. 
In consistency with the discussion of Fig. 23, the phase at which the supercurrent has a maximum, 4> ma x does not 
depend on the N-layer thickness, and it shifts to the right with increasing ds- 



VII. INTERFACE POTENTIALS 



Stimulated by recent work on the influence of interface barriers in SNS systems 14 , we studied this in more detail 
and for SN multilayers as well. Interface barriers can come out in practice as an effect of localized disorder at the 
interface or as a typical oxide layer in a point contact. 

A simple model of a (5-function potential at the interfaces introduced by Blonder et al. 13 , can be implemented in 
our formalism easily. The corresponding Hamiltonian for interfaces at positions Xj reads as 



'dx 2 



klx + ^WjSix-x,), 



(46) 



where Wj is the strength of the barrier and can be estimated using the transmission coefficient of the barrier 



T = 



1 



1 + mW 2 /Ti z E 



(47) 



In the presence of a ^-function barrier, the wave function is still given by equation (4), but the boundary conditions 
for the Green's function read now 



E 



SujGvju' j(Xj, 



x') = 0, 



(48) 



where 
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(49) 



We choose Wj = W the same for each SN interface. 

Before applying this to SN structures, we first look at the bound states of a SNS system at different choices of the 
barrier strength W. In Fig. 31 we show the LDOS of a SNS system characterized by dN = 10000 Bohr and L t = 13 
Bohr. In the absence of the interface potential, the SNS system LDOS has four bound states for E < A, corresponding 
to the modes (1,1), (1,2), (2,1), and (2,2). Due to the fact that L y = L z = L t , the states corresponding to the (1,2) 
and (2,1) modes are degenerate. Apparently, as we can notice, the presence of the 5-function potential favours the 
appearance of new bound states, due to the scattering with the interface potential. This implies that in the presence 
of a scattering potential even in AA the bound states split up. 

In the case of the multilayer, for which results are shown in Fig. 32, the deviation from the zero potential case is 
even more pronounced, as the quasi-periodicity of the dispersion relations (36) is perturbed by the interface barrier. 
As for the SNS system, new bound states appear and complicate the picture seen in Fig. 14, which describes the 
same multilayer, but in the absence of an interface barrier. 

In the limit of a large barrier strength W (W > 1 RydBohr), the S-layers decouple, so that the density of states for 
a multilayer becomes similar to the one of a SNS system. This can be seen in Fig. 33, in which we show the LDOS for 
a SN multilayer with ds = rfjv = 10000 and L t = 13, and for a SNS system with djy = 10000. The barrier strength 
is W — 10 RydBohr. Compared to the strength of a S layer, which is Ads = 1 RydBohr, this interface barrier is 10 
times larger. At such a large strength of the barrier, in the N-layer there are just bound states. To make it more clear, 
in Fig. 34 we show the contribution from each transverse mode to the LDOS of a SNS system, for energies up to 10 



10 



times the gap. The peaks corresponding to the same mode (n y ,n z ) have the same height and obey the dispersion 
relation for a 3-dimensional box, 




E + , = K + K + K=\^-^\ + m + ^ , (50) 



with 2n x + 1 > and = - (^) 2 - (^f ) 2 . Thus, for the (1,1) mode, the first peak has n x = 985 and 

is situated at the energy E — 2A, while the second peak has n x — 986 and occurs at E = 9.7 A. Similarly, using Eq. 
(50) for the modes (1,2) and (2,1), we obtain peaks for n x — 726 at E = 3.6A and for n x — 727 at E — 9.3A. For the 
mode (2,2) we get peaks for n x = 288, 289, 290, 291, and 292, at the energies E = 0.6A, 2.9A, 5.1A, 7.4A, and 9.7A 
respectively. 

The fact that the S-layers decouple in the limit of large barrier strength has also consequences on the phase 
dependence of the LDOS. We first show in Fig. 35 the LDOS of an SN multilayer without interface barrier for <f> = 
and cj) = tt. In this figure the solid curves of the Figs. 14 and 15 are shown in one picture. Clearly, the features of the 
LDOS, already discussed in Section III, are different for the two values of the phase <f). However, in the presence of 
an interface barrier, the picture changes. In Figs. 36 and 37 we show the LDOS for W = 1 and W = 10 respectively. 
When W = 1, the bound states occur almost at the same energies for both phases, and for W — 10 the LDOS almost 
coincide, as a result of a complete decoupling of the successive S-layers. This leads to a total suppression of the 
supercurrent at large values of W . 

The calculated supercurrent / for a SN multilayer at different barrier strengths W is shown in Fig. 38. Clearly, 
an interface barrier diminishes the supercurrent. The stronger the barrier, the smaller is the transmission probability 
through the interface. For values of W larger that 2 RydBohr the supercurrent is completely suppressed. A similar 
result was obtained in the recent study mentioned above, of the Josephson current in the much simpler SNS system 
having several insulating barriers 14 . 



VIII. CONCLUSIONS 



In this paper we discussed SN multilayer structures. In particular we show results for periodic infinite multilayers, 
represented by a Kronig-Pcnney superlattice model. By applying a Green's function formalism, we focussed first 
on the Andreev bound states and we studied the limitations of the Andreev approximation in relation to the finite 
transverse size of the systems. 

Further, we calculated the supercurrent through such a periodic SN multilayer. We completed this study using 
a selfconsistently calculated gap. Finally, including a <5-function potential at the interface, we derived results which 
account for possible barrier scattering at the interfaces. 

The results presented in this paper are meant to understand the physics which is behind SN multilayer structures. 
For our purpose it is more appropriate to investigate systems of very small transverse size because in such systems 
the effects of the breakdown of the Andreev approximation come out most clearly. However, at the moment there 
are no experimental data to which we can compare. For larger systems, more accessible to experiments, the physics 
remains the same, but their complexity could obscure some of the fundamental aspects we are looking at. 

Applications to intrinsic Josephson junctions 15 made from high-T c superconducting materials would require an 
extension of the present theory to the case of d-wave symmetry of the order parameter. 



APPENDIX: THE MATRICES A, B, AND C 



In this Appendix we show the structure of the equations (28) to (30) more explicitly The matrices A and C are 
highly singular and due to this, the solutions of Eqs. (28) and (29) are sparse matrices, which can be written 
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and 



X+ = 



( 6 
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Al4 
X + 
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-24 
X + 
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-34 

x + 
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-44 
X + 
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ooo 
6 6 6 
6 6 6 
6 6 6 



xf xf 
xf xf 
xf xf 
xf xf 
x 5 + 7 xf 
xf xf 

X 7 , 7 xf 

xf xf 



(A2) 



Substituting these matrices into Eqs. (28) and (29), we can reduce the set to solving two quadratic matrix equations 
for the 2x2 complex matrices x 3 ? and xf. 



1,33 



(A3) 



and 



X 77 Xf 

xf xf 



(A4) 



This appears to be equivalent to solving a system of 8 simultaneous equations with real coefficients. Mathematically, 
one can never predict the number of solutions. We solve this system numerically, by applying Newton's method, which 
requires an initial guess of the solution. Since we know that the solution is close to the Andreev approximation, wc 
give as an initial guess a diagonal matrix, namely the unitary matrix. This leads us to the physical solutions for X_ 
and X + which, by using equation (30), allows the calculation of the To matrix, which we need for the T 3 y matrix 
(27) and the Green's function (31). Finally, we are able to make use of the expressions (2) and (3) and calculate the 
LDOS and the supercurrent of a periodic SN multilayer. 

In Andreev approximation some matrix elements are zero and further simplifications can be made. By neglecting 
the ordinary reflections of the quasiparticles at the SN interfaces, some of the t-matrices are equal to zero. This results 
in a more simple form for the matrices A and C 



A = 



\ 

fl 2 4 



/ 



(A5) 



where 



a 2 4 = 



-ik , . 04 

+3 ^ 



f-++ 







;+-+- 



j— ik . .a_i_ 

d +j e +, 

ik, aa 



- j— ik . .t 

d +j e +* 



(A6) 



where 



Cis 

c = i S S S ;: ! ■ 





. _ / i+t- + dt 3 e- lk+ -i a - t+Z'+dz/^ 

013 ~ ^ ii+- + dt je ik >- iiz- + dz je - ik: 



+ " 3 ■ (A8) 

t+Z- + dZ J e- lk -i a - J 
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By consequence, the solutions x 3 ? and x A } get then a diagonal form. 



and 



*t={^ (A10) 

This allows us to decouple the equations for X^p and X^p into two quadratic equations which can be solved directly. 
The same holds for x+ . Combining the two possible solutions for X^_ 5 and X 66 with the two possible solutions for 

X^J and X^ 8 , one obtains four mathematical solutions for X and X+. Two of the solutions are complementary and 
lead to a zero value for the LDOS. Using the other two solutions, one gets either the positive physical value for the 
LDOS, or the same value with the opposite sign. We use this criterion to distinguish the physical solution from the 
four mathematically possible solutions. 
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FIG. 1. Illustration of a Superconductor/Normal metal multilayer with finite transverse widths. 
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FIG. 2. The Kronig-Penney model for the pair potential, used in Ref.* 
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FIG. 3. Schematic reDresentation of the processes accounted for bv the matrices Aj , and Cj 
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FIG. 4. LDOS for a SN multilayer (ds = 50000 Bohr and djv = 1000 Bohr) and for a bar shaped superconductor (dashed 
line) . 
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FIG. 5. LDOS for a SN multilayer (ds = 50000 Bohr and d N = 1000 Bohr) in the N layer (solid line) and S layer (dashed 
line) . 



LDOS inside the N layer 
LDOS inside the S layer 



d s =50000 Bohr 
d N =2000 Bohr 




FIG. 6. LDOS for a SN multilayer (ds = 50000 Bohr and d N = 2000 Bohr) in the N layer (solid line) and S layer (dashed 
line) . 




FIG. 7. LDOS for a SN multilayer (ds = 50000 Bohr and d N = 4000 Bohr) in the N layer (solid line) and S layer (dashed 
line) . 
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FIG. 8. LDOS for a SN multilayer (d s = 50000 Bohr and d N = 10000 Bohr) in the N layer (solid line) and S layer (dashed 
line) . 



SN multilayer, in the N layer 
SNS system, in the N layer 



d s =50000 Bohr 
d N =100()0 Bohr 




FIG. 9. LDOS for a SN multilayer (d s = 50000 Bohr and d N = 10000 Bohr) and for a SNS system (dashed line), calculated 
inside the N layer. 
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FIG. 10. LDOS for a SN multilayer (ds = 30000 Bohr and d N = 10000 Bohr) and for a SNS system (dashed line), calculated 
inside the N layer. 
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FIG. 11. LDOS for a SN multilayer (d s = 10000 Bohr and d N = 10000 Bohr) and for a SNS system (dashed line), calculated 
inside the N layer. 




FIG. 12. Contributions to the LDOS of a SN multilayer (ds = djv = 10000 Bohr) in the middle of the N layer from the mode 
(1,2), for two choices of the phase of the pair potential, <f> = and (j> — ir. 
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FIG. 13. Contributions to the LDOS of a SN multilayer (ds = djv = 10000 Bohr) in the middle of the N layer from the 
modes (1,1), (1,2), (2,1), and (2,2). The phase of the pair potential is <f> = 0. 
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E/A 

FIG. 14. LDOS for a SN multilayer (ds = djv = 10000 Bohr) in the middle of the N layer (solid line) and S layer (dashed 
line). The phase of the pair potential is d> — 0. 




FIG. 15. LDOS for a SN multilayer (ds — djv = 10000 Bohr) in the N layer (solid line) and S layer (dashed line). The phase 
of the pair potential is <f> — n. 




FIG. 16. LDOS calculated in the N and S layers of a SN multilayer, at a higher transverse size, Lt = 130 Bohr. 
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FIG. 17. LDOS calculated in a bar shaDcd S material, at E = 5A. at different transverse widths. 




FIG. 18. LDOS calculated in a SN multilayer, for a width L t = 12.566371 Bohr close to a critical width Lf , at different posi- 
tions x with respect to the S/N interface, inside the N layer. d N = d$ = 10000 Bohr. For comparison, Andreev approximation 
is represented with dotted line. 
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FIG. 19. LDOS calculated in both S and N layers of a SN multilayer, with djv = ds = 4000 Bohr, as well as in a SNS system, 
at a transverse width L t = 12.5676 Bohr. 
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FIG. 20. LDOS calculated in Andreev approximation in the N layer of a SN multilayer, with djv = ds = 4000 Bohr, as well 
as in a SNS system, at a transverse width L f — 12.5676 Bohr. 




FIG. 21. The dependence of the supercurrent I on the phase of the pair potential <f>. d,N = ds = 10000 Bohr and the 
normalization factor 7o = eA/h. 
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FIG. 22. The dependence of the supercurrent I on the transverse length L t . Io = eA/h. 
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FIG. 23. The dependence of the supercurrent / on djv and ds for L t — 13 Bohr. Io = eA/h. 
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FIG. 24. Selfconsistent gap calculations for systems with transverse width Lt=30 Bohr. 
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FIG. 25. Selfconsistent gap calculations for systems with transverse width L t = 100 Bohr. 
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FIG. 26. Selfconsistent gap calculations for homogeneous S bars of different thicknesses L t , at T = 0, for two different cut-off 
methods. 
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FIG. 27. Selfconsistent gap calculations for an homogeneous S bar of L t =30 Bohr, using momenta cut-off. 
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FIG. 29. Selfconsistent calculation of the supercurrent as a function of temperature for a SN multilayer with ds = d]v=10000 
and Lt=30. J = eA s /h and A s = 0.8 x 1(T 4 . 
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FIG. 30. Selfconsistent calculation of the supercurrent as a function of phase for SN multilayers with d$ = djv = 10000 
and ds = djv=4000, and corresponding SNS systems with djv=10000 and djv=4000 respectively. L t —30, Io = eA s /h, and 
A s = 0.8 x 10~ 4 . 




FIG. 31. LDOS for a SNS svstem at different barrier streneths W. = 10000, L t = 13. 
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FIG. 32. LDOS for a SN multilayer at different barrier strengths W. d s = d N = 10000, L t = 13. 
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FIG. 33. LDOS in the N-layer of a SN multilayer and of a SNS system, at W = 10 RydBohr. d s = d N = 10000, L t = 13. 
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FIG. 34. The contributions of the transverse modes to the LDOS of a SNS system with W = 10. d N = 10000, L t = 13. 




FIG. 35. LDOS in the N-layer of a SN multilayer at (j> = and 7r, and without interface barrier, ds = d,N = 10000, L t = 13. 
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FIG. 36. LDOS in the N-layer of a SN multilayer at <j> = and tt, and with W = 1. ds = d N = 10000, L t = 13. 
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FIG. 37. LDOS in the N-layer of a SN multilayer at <t> = and tt, and with W = 10. d s = d N = 10000, L t = 13. 




FIG. 38. Supercurrent for a SN multilayer at different potentials, djv = 10000, L t = 13, and Iq — eA/h. 
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